How does urbanization influence ibis behaviors relevant to transmission?

Specifically, how does variation in flock size and frequency of provisioning influence these behaviors?

Data:

10-minute focal follow observations of birds at 5 urban parks in South Florida. Follows conducted during summer 2019.

Focal follows were only conducted in urban parks

Note: Focal follows were only conducted when people were not actively feeding the birds so foraging means probing on the ground for food. The few instances of consuming human provided food are cases where ibis were foraging ‘naturally’ and happened upon scraps of anthropogenic food that was already on the ground.

Sample Sizes / Data Structure

Observations consist of 168 focal follows lasting > 2 minutes. Five urban parks were samples for ~30 focals per site. Sites were visited 2-4 times each.

Provisioning estimates: Each visit, # of people and groups feeding birds were recorded for duration of visit. Provisioning was estimated by counting the number of groups that fed the birds per visit and dividing by total time I was at the site to get # provisioning events/ hour. Each focal follow was paired with the average # of provisioning events per hour on day/site of the focal follow.

Flock size estimates: Every ~30 minutes, flock size was recorded. I paired each focal follow with the nearest estimate of flock size.

Birds seem to spend the most time foraging, grooming, and being vigilant.

I am going to re-group these variables according to their relevance to different modes of parasite transmission (see different sections below).

Summary Statistics

How many focal follows lasted > 2 minutes? *important- when we add in predictors later, we lose data for 3 follows because they lack corresponding predictor data.

## [1] 150

What was the total amount of time spent observing ibis? (this only includes obs >2 minutes aka the ones used in the analysis)

Total obs time in seconds:

## [1] 79994

Total obs time in minutes:

## [1] 1333.233

Mean obs time:

## [1] 533.2933

How many different behaviors were observed on average for each ibis during the observation?

##    focalID           n_behaviors   
##  Length:150         Min.   :2.000  
##  Class :character   1st Qu.:3.000  
##  Mode  :character   Median :4.000  
##                     Mean   :3.853  
##                     3rd Qu.:5.000  
##                     Max.   :6.000

How many birds showed 2 or more behaviors?

nrow(sum3.1[sum3.1$n_behaviors >= 2,]) #all birds observed for at least 2 minutes did 2 or more behaviors 
## [1] 150
nrow(sum3.1[sum3.1$n_behaviors < 2,])
## [1] 0

Most commonly observed behaviors observed?

How many birds foraged?

sum4 <- simp2

sum4f <- subset(sum4, Activity=="Foraging" & activity_count > 0)
summary(sum4f)
##    focalID          focalduration     Activity         activity_count 
##  Length:110         Min.   :139.0   Length:110         Min.   :  4.0  
##  Class :character   1st Qu.:592.0   Class :character   1st Qu.: 53.5  
##  Mode  :character   Median :596.0   Mode  :character   Median :148.0  
##                     Mean   :537.1                      Mean   :212.6  
##                     3rd Qu.:597.0                      3rd Qu.:362.0  
##                     Max.   :598.0                      Max.   :585.0

How many birds groomed?

sum4g <- subset(sum4, Activity=="Grooming" & activity_count > 0)
summary(sum4g)
##    focalID          focalduration     Activity         activity_count 
##  Length:113         Min.   :122.0   Length:113         Min.   :  2.0  
##  Class :character   1st Qu.:557.0   Class :character   1st Qu.: 22.0  
##  Mode  :character   Median :596.0   Mode  :character   Median :140.0  
##                     Mean   :536.5                      Mean   :170.8  
##                     3rd Qu.:597.0                      3rd Qu.:275.0  
##                     Max.   :598.0                      Max.   :572.0

How many birds were vigilant?

sum4v <- subset(sum4, Activity=="Vigilant" & activity_count > 0)
summary(sum4v)
##    focalID          focalduration     Activity         activity_count  
##  Length:140         Min.   :122.0   Length:140         Min.   :  4.00  
##  Class :character   1st Qu.:579.8   Class :character   1st Qu.: 41.75  
##  Mode  :character   Median :596.0   Mode  :character   Median :132.50  
##                     Mean   :539.2                      Mean   :178.65  
##                     3rd Qu.:597.0                      3rd Qu.:269.00  
##                     Max.   :598.0                      Max.   :591.00

What behaviors did ibis spend the most time doing?

sum5<-simp2

sum5 <- sum5 %>% group_by(Activity) %>%dplyr::summarise(sum(activity_count))
sum5
## # A tibble: 7 × 2
##   Activity `sum(activity_count)`
##   <chr>                    <dbl>
## 1 Drinking                   243
## 2 Foraging                 23387
## 3 Grooming                 19305
## 4 Other                      812
## 5 Resting                   2934
## 6 Vigilant                 25011
## 7 Walking                   8302

Which behaviors were observed most frequently in ibis?

sum6<-simp2

sum6<- as.data.frame(subset(simp2, activity_count > 0))

sum6.1 <- sum6 %>% group_by(Activity) %>% 
 dplyr::summarise(n_focals=n(),
            .groups = 'drop')

sum6.1
## # A tibble: 7 × 2
##   Activity n_focals
##   <chr>       <int>
## 1 Drinking       19
## 2 Foraging      110
## 3 Grooming      113
## 4 Other          61
## 5 Resting        20
## 6 Vigilant      140
## 7 Walking       115

Visualizations of predictors

Predictors:

What is the distribution of provisioning frequency? What about average flock size? Do we see any relationship between # feeding events and average flock size?

Note: Here I’m showing avg flock size and nearest flock size

# Provisioning and Flock Density

ggplot(predictors, aes(x=feedingperhour, y=flocksize_nearest)) +
  geom_point() 
## Warning: Removed 3 rows containing missing values (`geom_point()`).

ggplot(predictors, aes(x=feedingperhour, y=log(flocksize_nearest))) +
  geom_point() 
## Warning: Removed 3 rows containing missing values (`geom_point()`).

lm1 <- lm(predictors$flocksize_nearest~predictors$feedingperhour)
summary(lm1)
## 
## Call:
## lm(formula = predictors$flocksize_nearest ~ predictors$feedingperhour)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -46.368 -31.945 -16.797   5.128 169.065 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 40.111      7.049    5.69 6.58e-08 ***
## predictors$feedingperhour    9.647      8.178    1.18     0.24    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 53.47 on 148 degrees of freedom
##   (3 observations deleted due to missingness)
## Multiple R-squared:  0.009315,   Adjusted R-squared:  0.002621 
## F-statistic: 1.392 on 1 and 148 DF,  p-value: 0.24

Model visualizations and Model Selection

Standardize predictor variables

#library(mosaic)  # standardizing variables

ForagePreen_wide <- 
  ForagePreen_wide %>% 
  mutate(avg_flock_size_s = scale(avg_flock_size))%>% 
  mutate(flocksize_nearest_s = scale(flocksize_nearest))%>%
  mutate(feedingperhour_s = scale(feedingperhour))

Vig_wide <- 
  Vig_wide %>% 
  mutate(avg_flock_size_s = scale(avg_flock_size))%>% 
  mutate(flocksize_nearest_s = scale(flocksize_nearest))%>%
  mutate(feedingperhour_s = scale(feedingperhour))

Foraging models:

Main Q:

Does flock size, provisioning, both, or neither influence foraging time?

Foraging includes foraging on soil, in water, on human objects, and consuming human provided food.

Note: Focal follows were only conducted when people were not actively feeding the birds so foraging means probing on the ground for food. The few instances of consuming human provided food are cases where ibis were foraging ‘naturally’ and happened upon scraps of anthropogenic food that was already on the ground.

Exploratory Visualizations

Model selection:

Count data: Does poisson or negative binomial fit better? Build global model (two predictors + interactions) and compare AIC.

pm1 = glmmTMB(Foraging~flocksize_nearest_s*feedingperhour_s + offset(log(duration)), ForagePreen_wide, family=poisson)

nbm1 = glmmTMB(Foraging~flocksize_nearest_s*feedingperhour_s  +offset(log(duration)), ForagePreen_wide, family=nbinom2)

zipm1 = glmmTMB(Foraging~flocksize_nearest_s*feedingperhour_s  +offset(log(duration)), zi=~flocksize_nearest_s*feedingperhour_s, ForagePreen_wide, family=poisson)

zipm1.int = glmmTMB(Foraging~flocksize_nearest_s*feedingperhour_s  +offset(log(duration)), zi=~1, ForagePreen_wide, family=poisson)

zinbm1 = glmmTMB(Foraging~flocksize_nearest_s*feedingperhour_s +offset(log(duration)), zi=~flocksize_nearest_s*feedingperhour_s, ForagePreen_wide, family=nbinom2)

zinbm1.int = glmmTMB(Foraging~flocksize_nearest_s*feedingperhour_s +offset(log(duration)), zi=~1, ForagePreen_wide, family=nbinom2) #intercept only zero inflated term (equally applies zero inflated to whole model)


AICtab(pm1,nbm1,zipm1,zipm1.int,zinbm1,zinbm1.int) #zinbm1.int is best fit 
##            dAIC    df
## zinbm1.int     0.0 6 
## zinbm1         2.5 9 
## nbm1          59.9 5 
## zipm1.int  14290.9 5 
## zipm1      14293.2 8 
## pm1        27285.8 4

Is this the case for other candidate models (not just the global model)?

#just flock size
pm2 = glmmTMB(Foraging~flocksize_nearest_s + offset(log(duration)), ForagePreen_wide, family=poisson)

nbm2 = glmmTMB(Foraging~flocksize_nearest_s + offset(log(duration)), ForagePreen_wide, family=nbinom2)

zipm2 = glmmTMB(Foraging~flocksize_nearest_s + offset(log(duration)), zi=~flocksize_nearest_s, ForagePreen_wide, family=poisson)

zinbm2 = glmmTMB(Foraging~flocksize_nearest_s + offset(log(duration)), zi=~flocksize_nearest_s, ForagePreen_wide, family=nbinom2)

zinbm2.int = glmmTMB(Foraging~flocksize_nearest_s + offset(log(duration)), zi=~1, ForagePreen_wide, family=nbinom2)

AICtab(pm2,nbm2,zipm2,zinbm2,zinbm2.int) #the two zinbm models are the best 
##            dAIC    df
## zinbm2.int     0.0 4 
## zinbm2         1.8 5 
## nbm2          58.6 3 
## zipm2      15518.7 4 
## pm2        29415.6 2
#just provisioning
pm3 = glmmTMB(Foraging~feedingperhour_s + offset(log(duration)), ForagePreen_wide, family=poisson)

nbm3 = glmmTMB(Foraging~feedingperhour_s + offset(log(duration)), ForagePreen_wide, family=nbinom2)

zipm3 = glmmTMB(Foraging~feedingperhour_s + offset(log(duration)), zi=~feedingperhour_s, ForagePreen_wide, family=poisson)

zinbm3 = glmmTMB(Foraging~feedingperhour_s + offset(log(duration)), zi=~feedingperhour_s, ForagePreen_wide, family=nbinom2)

zinbm3.int = glmmTMB(Foraging~feedingperhour_s + offset(log(duration)), zi=~1, ForagePreen_wide, family=nbinom2)

AICtab(pm3,nbm3,zipm3,zinbm3,zinbm3.int) #the two zinbm models are the best
##            dAIC    df
## zinbm3.int     0.0 4 
## zinbm3         0.2 5 
## nbm3          57.8 3 
## zipm3      14817.4 4 
## pm3        28058.0 2
#Flock size and provisioning 
pm4 = glmmTMB(Foraging~flocksize_nearest_s+feedingperhour_s + offset(log(duration)), ForagePreen_wide, family=poisson)

nbm4 = glmmTMB(Foraging~flocksize_nearest_s+feedingperhour_s  +offset(log(duration)), ForagePreen_wide, family=nbinom2)

zipm4 = glmmTMB(Foraging~flocksize_nearest_s+feedingperhour_s  +offset(log(duration)), zi=~flocksize_nearest_s+feedingperhour_s, ForagePreen_wide, family=poisson)

zinbm4 = glmmTMB(Foraging~flocksize_nearest_s+feedingperhour_s +offset(log(duration)), zi=~flocksize_nearest_s+feedingperhour_s, ForagePreen_wide, family=nbinom2)

zinbm4.int = glmmTMB(Foraging~flocksize_nearest_s+feedingperhour_s +offset(log(duration)), zi=~1, ForagePreen_wide, family=nbinom2) #intercept only zero inflated term

AICtab(pm4,nbm4,zipm4,zinbm4,zinbm4.int) #the two zinbm models are the best 
##            dAIC    df
## zinbm4.int     0.0 5 
## zinbm4         2.0 7 
## nbm4          60.2 4 
## zipm4      14327.2 6 
## pm4        27600.9 3

Yes- the zero inflated models always outperform the others.

I can also double check dispersion and zero inflation in my models using Dharma:

#check fit of poisson model: 
simulationOutput.pm1 <- simulateResiduals(fittedModel = pm1)
plot(simulationOutput.pm1)

#overdispersed or underdispersed?
testDispersion(simulationOutput.pm1,alternative = "greater") #test for overdispersion 

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 186.35, p-value < 2.2e-16
## alternative hypothesis: greater
testDispersion(simulationOutput.pm1,alternative = "less") #test for underdispersion 

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 186.35, p-value = 1
## alternative hypothesis: less
#poisson model is overdispersed 

#zero inflated?
testZeroInflation(simulationOutput.pm1) #yes zero inflated 

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = Inf, p-value < 2.2e-16
## alternative hypothesis: two.sided
#Poisson model is overdispersed and zero inflated- will the neg binom model do better?
simulationOutput.nbm <- simulateResiduals(fittedModel = nbm1)
plot(simulationOutput.nbm)

#overdispersed or underdispersed?
testDispersion(simulationOutput.nbm,alternative = "greater") #test for overdispersion 

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.29045, p-value = 1
## alternative hypothesis: greater
testDispersion(simulationOutput.nbm,alternative = "less") #test for underdispersion 

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.29045, p-value < 2.2e-16
## alternative hypothesis: less
#This model is underdispersed

testZeroInflation(simulationOutput.nbm) #still slightly zero inflated 

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.5642, p-value = 0.016
## alternative hypothesis: two.sided
#Now test a zero inflated neg binom model (look at zinbm1 and zinbm1.int)
simulationOutput.zinbm1 <- simulateResiduals(fittedModel = zinbm1)
plot(simulationOutput.zinbm1)

simulationOutput.zinbm1.int <- simulateResiduals(fittedModel = zinbm1.int)
plot(simulationOutput.zinbm1.int)

testDispersion(simulationOutput.zinbm1,alternative = "greater") #test for overdispersion 

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.6559, p-value = 0.924
## alternative hypothesis: greater
testDispersion(simulationOutput.zinbm1,alternative = "less") #test for underdispersion 

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.6559, p-value = 0.076
## alternative hypothesis: less
testDispersion(simulationOutput.zinbm1.int,alternative = "greater") #test for overdispersion 

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.67235, p-value = 0.9
## alternative hypothesis: greater
testDispersion(simulationOutput.zinbm1.int,alternative = "less") #test for underdispersion 

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.67235, p-value = 0.1
## alternative hypothesis: less
#disperson is fine for both of these, what about zero-inflation?
testZeroInflation(simulationOutput.zinbm1)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0068, p-value = 1
## alternative hypothesis: two.sided
testZeroInflation(simulationOutput.zinbm1.int)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0021, p-value = 1
## alternative hypothesis: two.sided
#Both of these models are good- no dispersion or zero inflation issues

These results suggest that the zero inflated models make the most sense (they fix the dispersion and zero inflation issues). ZI models with zi=1 and zi= conditional fixed effects are both sufficient in fixing these issues.

For these data, I would think most zeros are false zeros that are due to experimental design AKA 0 foraging does not mean the bird never forages but rather it was not observed for long enough to observe foraging. I don’t think these false zeros could be explained by my predictors. Because of this, I think I should apply a single zero inflation parameter to my entire model (ziformula~1).

Now I will move forward with zero inflated negative binomial models with a single zero inflated term applied to all models.

First build models without any random effects:

#intercept only
fmod0.3 <-glmmTMB(Foraging~1+ offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=ForagePreen_wide)
summary(fmod0.3)
##  Family: nbinom2  ( log )
## Formula:          Foraging ~ 1 + offset(log(duration))
## Zero inflation:            ~1
## Data: ForagePreen_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1571.0   1580.1   -782.5   1565.0      147 
## 
## 
## Dispersion parameter for nbinom2 family (): 0.987 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.91841    0.09627   -9.54   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.0340     0.1881  -5.497 3.85e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#provisioning
fmod11 <-glmmTMB(Foraging~feedingperhour_s+offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=ForagePreen_wide)
summary(fmod11)
##  Family: nbinom2  ( log )
## Formula:          Foraging ~ feedingperhour_s + offset(log(duration))
## Zero inflation:            ~1
## Data: ForagePreen_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1567.9   1579.9   -779.9   1559.9      146 
## 
## 
## Dispersion parameter for nbinom2 family (): 1.03 
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -0.95398    0.09454 -10.091   <2e-16 ***
## feedingperhour_s -0.19827    0.08581  -2.311   0.0209 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.0324     0.1878  -5.497 3.86e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#flock size
fmod12 <-glmmTMB(Foraging~flocksize_nearest_s+offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=ForagePreen_wide)
summary(fmod12)
##  Family: nbinom2  ( log )
## Formula:          Foraging ~ flocksize_nearest_s + offset(log(duration))
## Zero inflation:            ~1
## Data: ForagePreen_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1567.6   1579.6   -779.8   1559.6      146 
## 
## 
## Dispersion parameter for nbinom2 family (): 1.03 
## 
## Conditional model:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -0.94968    0.09438  -10.06   <2e-16 ***
## flocksize_nearest_s -0.26320    0.10240   -2.57   0.0102 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.0343     0.1881    -5.5  3.8e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#provisioning + flock size
fmod13 <-glmmTMB(Foraging~flocksize_nearest_s+feedingperhour_s+offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=ForagePreen_wide)
summary(fmod13)
##  Family: nbinom2  ( log )
## Formula:          
## Foraging ~ flocksize_nearest_s + feedingperhour_s + offset(log(duration))
## Zero inflation:            ~1
## Data: ForagePreen_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1565.2   1580.3   -777.6   1555.2      145 
## 
## 
## Dispersion parameter for nbinom2 family (): 1.06 
## 
## Conditional model:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -0.97975    0.09296 -10.539   <2e-16 ***
## flocksize_nearest_s -0.24333    0.10301  -2.362   0.0182 *  
## feedingperhour_s    -0.18472    0.08670  -2.131   0.0331 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.0330     0.1878  -5.499 3.82e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fmod13.1 <-glmmTMB(Foraging~flocksize_nearest_s+feedingperhour_s+offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=ForagePreen_wide)
summary(fmod13.1)
##  Family: nbinom2  ( log )
## Formula:          
## Foraging ~ flocksize_nearest_s + feedingperhour_s + offset(log(duration))
## Zero inflation:            ~1
## Data: ForagePreen_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1565.2   1580.3   -777.6   1555.2      145 
## 
## 
## Dispersion parameter for nbinom2 family (): 1.06 
## 
## Conditional model:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -0.97975    0.09296 -10.539   <2e-16 ***
## flocksize_nearest_s -0.24333    0.10301  -2.362   0.0182 *  
## feedingperhour_s    -0.18472    0.08670  -2.131   0.0331 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.0330     0.1878  -5.499 3.82e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#provisioning x flock size 
fmod14 <-glmmTMB(Foraging~flocksize_nearest_s*feedingperhour_s+offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=ForagePreen_wide)
summary(fmod14)
##  Family: nbinom2  ( log )
## Formula:          
## Foraging ~ flocksize_nearest_s * feedingperhour_s + offset(log(duration))
## Zero inflation:            ~1
## Data: ForagePreen_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1567.1   1585.2   -777.6   1555.1      144 
## 
## 
## Dispersion parameter for nbinom2 family (): 1.07 
## 
## Conditional model:
##                                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                          -0.97711    0.09345 -10.456   <2e-16 ***
## flocksize_nearest_s                  -0.24949    0.10501  -2.376   0.0175 *  
## feedingperhour_s                     -0.20032    0.09952  -2.013   0.0441 *  
## flocksize_nearest_s:feedingperhour_s -0.05511    0.17401  -0.317   0.7515    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.0333     0.1879  -5.499 3.81e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AICtab(fmod0.3,fmod11,fmod12,fmod13,fmod14)
##         dAIC df
## fmod13  0.0  5 
## fmod14  1.9  6 
## fmod12  2.4  4 
## fmod11  2.7  4 
## fmod0.3 5.8  3

fmod13 (flock size and feeding) and fmod14 (flock size * feeding) are best fitting models.

Now build the same models but include site as a random effect and see if same models are best fitting. I think site makes sense biologically but it only has 4 levels so questionable as a RE.

#intercept only
fmod0.2 <-glmmTMB(Foraging~1+ (1|site) + offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=ForagePreen_wide)
summary(fmod0.2)
##  Family: nbinom2  ( log )
## Formula:          Foraging ~ 1 + (1 | site) + offset(log(duration))
## Zero inflation:            ~1
## Data: ForagePreen_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1570.3   1582.4   -781.2   1562.3      146 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.08102  0.2846  
## Number of obs: 150, groups:  site, 5
## 
## Dispersion parameter for nbinom2 family (): 1.05 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.9674     0.1599  -6.052 1.43e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -1.033      0.188  -5.498 3.84e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#provisioning
fmod15 <-glmmTMB(Foraging~feedingperhour_s+(1|site)+offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=ForagePreen_wide)
summary(fmod15)
##  Family: nbinom2  ( log )
## Formula:          
## Foraging ~ feedingperhour_s + (1 | site) + offset(log(duration))
## Zero inflation:            ~1
## Data: ForagePreen_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1569.8   1584.9   -779.9   1559.8      145 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.01202  0.1097  
## Number of obs: 150, groups:  site, 5
## 
## Dispersion parameter for nbinom2 family (): 1.04 
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -0.9640     0.1126  -8.561   <2e-16 ***
## feedingperhour_s  -0.1861     0.1062  -1.752   0.0798 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.0327     0.1879  -5.497 3.86e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#flock size
fmod16 <-glmmTMB(Foraging~flocksize_nearest_s+(1|site)+offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=ForagePreen_wide)
summary(fmod16)
##  Family: nbinom2  ( log )
## Formula:          
## Foraging ~ flocksize_nearest_s + (1 | site) + offset(log(duration))
## Zero inflation:            ~1
## Data: ForagePreen_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1569.6   1584.6   -779.8   1559.6      145 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev. 
##  site   (Intercept) 4.094e-07 0.0006398
## Number of obs: 150, groups:  site, 5
## 
## Dispersion parameter for nbinom2 family (): 1.03 
## 
## Conditional model:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -0.94967    0.09439 -10.062   <2e-16 ***
## flocksize_nearest_s -0.26319    0.10244  -2.569   0.0102 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.0343     0.1881    -5.5  3.8e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#provisioning + flock size
fmod17 <-glmmTMB(Foraging~flocksize_nearest_s+feedingperhour_s+(1|site)+offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=ForagePreen_wide)
summary(fmod17)
##  Family: nbinom2  ( log )
## Formula:          
## Foraging ~ flocksize_nearest_s + feedingperhour_s + (1 | site) +  
##     offset(log(duration))
## Zero inflation:            ~1
## Data: ForagePreen_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1567.2   1585.3   -777.6   1555.2      144 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev. 
##  site   (Intercept) 2.283e-09 4.779e-05
## Number of obs: 150, groups:  site, 5
## 
## Dispersion parameter for nbinom2 family (): 1.06 
## 
## Conditional model:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -0.97975    0.09296 -10.539   <2e-16 ***
## flocksize_nearest_s -0.24333    0.10301  -2.362   0.0182 *  
## feedingperhour_s    -0.18472    0.08670  -2.131   0.0331 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.0330     0.1878  -5.499 3.82e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#provisioning x flock size 
fmod18 <-glmmTMB(Foraging~flocksize_nearest_s*feedingperhour_s+(1|site)+offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=ForagePreen_wide)
summary(fmod18)
##  Family: nbinom2  ( log )
## Formula:          
## Foraging ~ flocksize_nearest_s * feedingperhour_s + (1 | site) +  
##     offset(log(duration))
## Zero inflation:            ~1
## Data: ForagePreen_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1569.1   1590.2   -777.6   1555.1      143 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev. 
##  site   (Intercept) 2.117e-09 4.601e-05
## Number of obs: 150, groups:  site, 5
## 
## Dispersion parameter for nbinom2 family (): 1.07 
## 
## Conditional model:
##                                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                          -0.97710    0.09345 -10.455   <2e-16 ***
## flocksize_nearest_s                  -0.24949    0.10501  -2.376   0.0175 *  
## feedingperhour_s                     -0.20032    0.09952  -2.013   0.0441 *  
## flocksize_nearest_s:feedingperhour_s -0.05511    0.17401  -0.317   0.7515    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.0333     0.1879  -5.499 3.81e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AICtab(fmod0.2,fmod15,fmod16,fmod17,fmod18) #fmod17 and fmod18 = best fit
##         dAIC df
## fmod17  0.0  6 
## fmod18  1.9  7 
## fmod16  2.4  5 
## fmod15  2.6  5 
## fmod0.2 3.1  4

Models with flock size and feeding or flock size * feeding are still the best fitting models (fmod17 and fmod18 here).

Summary estimates aren’t too different for the top performing models regardless of whether they include site as a random effect or not.

Model diagnostics

Plotting predicted values:

Plots made using ggeffects: https://strengejacke.github.io/ggeffects/index.html

The ggeffects package computes estimated marginal means (predicted values) for the response, at the margin of specific values or levels from certain model terms, i.e. it generates predictions by a model by holding the non-focal variables constant and varying the focal variable(s).

First try this for the models that don’t include random effects (fmod13 and fmod14)

#predict(fmod13,ForagePreen_wide,se.fit=TRUE) #couold use predict fxn but ggpredict makes it easier to graph:

#just one variable in model
fmod13_predict <- ggpredict(fmod13, terms = "feedingperhour_s")
ggplot(fmod13_predict, aes(x, predicted)) +
  geom_line() +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .1)

plot(ggpredict(fmod13, terms = "feedingperhour_s"))

fmod13_predict2 <- ggpredict(fmod13, terms = c("feedingperhour_s","flocksize_nearest_s")) #both model terms

ggplot(fmod13_predict2, aes(x = x, y = predicted,color=group)) +
  geom_line() 

plot(fmod13_predict2)

#with interaction:
#note the interaction is not significant in this model

#just one variable in model
fmod14_predict <- ggpredict(fmod14, terms = "feedingperhour_s")
ggplot(fmod14_predict, aes(x, predicted)) +
  geom_line() +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .1)

plot(ggpredict(fmod14, terms = "feedingperhour_s"))

fmod14_predict2 <- ggpredict(fmod14, terms = c("feedingperhour_s","flocksize_nearest_s")) #both model terms

ggplot(fmod14_predict2, aes(x = x, y = predicted,color=group)) +
  geom_line() 

plot(fmod14_predict2)

The models look very similar for both fmod13 and fmod14. When you look at the model summary for fmod14, the interaction is not significant. The interaction adds another parameter to the model (this interaction is a ‘pretending variable’). So they’re both supported but fmod14 has more parameters.

Now try this for models with site as random effect (fmod17 and fmod18)

#just one variable in model
fmod17_predict <- ggpredict(fmod17, terms = "feedingperhour_s")
ggplot(fmod17_predict, aes(x, predicted)) +
  geom_line() +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .1)

plot(ggpredict(fmod17, terms = "feedingperhour_s"))

fmod17_predict2 <- ggpredict(fmod17, terms = c("feedingperhour_s","flocksize_nearest_s")) #both model terms

ggplot(fmod17_predict2, aes(x = x, y = predicted,color=group)) +
  geom_line() 

plot(fmod17_predict2)

#with interaction:
#note the interaction is not significant in this model

#just one variable in model
fmod18_predict <- ggpredict(fmod18, terms = "feedingperhour_s")
ggplot(fmod14_predict, aes(x, predicted)) +
  geom_line() +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .1)

plot(ggpredict(fmod18, terms = "feedingperhour_s")) 

plot(ggpredict(fmod18, terms = "flocksize_nearest_s"))

fmod18_predict2 <- ggpredict(fmod18, terms = c("feedingperhour_s","flocksize_nearest_s")) #both model terms

ggplot(fmod18_predict2, aes(x = x, y = predicted,color=group)) +
  geom_line() 

plot(fmod18_predict2)

For the paper, showing model outputs for fmod17 (flock size + provisioning) might be useful. Choosing fmod17 because it doesn’t include the interaction in fmod18 which isn’t significant and does include site as random effect (as opposed to fmod13)

Make a 2 panel figure showing Foraging on Y axis and 1) provisioning and 2) flock size on X axis.

head(fmod17_predict) #data frame for provisioning
## # Predicted counts of Foraging
## 
## feedingperhour_s | Predicted |           95% CI
## -----------------------------------------------
##            -1.26 |    252.83 | [192.44, 332.17]
##            -0.75 |    229.93 | [185.30, 285.32]
##            -0.41 |    216.15 | [178.51, 261.71]
##            -0.33 |    212.78 | [176.44, 256.62]
##             0.49 |    182.75 | [148.85, 224.38]
##             0.54 |    181.09 | [146.90, 223.23]
## 
## Adjusted for:
## * flocksize_nearest_s =   0.00
## *            duration = 533.29
## *                site = NA (population-level)
fmod17_predict3 <- ggpredict(fmod17, terms = "flocksize_nearest_s")
head(fmod17_predict3)
## # Predicted counts of Foraging
## 
## flocksize_nearest_s | Predicted |           95% CI
## --------------------------------------------------
##               -1.00 |    255.35 | [195.34, 333.80]
##               -0.80 |    243.22 | [191.35, 309.15]
##               -0.60 |    231.67 | [186.71, 287.46]
##               -0.40 |    220.67 | [181.19, 268.74]
##               -0.20 |    210.18 | [174.61, 253.01]
##                0.00 |    200.20 | [166.85, 240.21]
## 
## Adjusted for:
## * feedingperhour_s =   0.00
## *         duration = 533.29
## *             site = NA (population-level)
#merge these two dfs and then facet
#or just make each one pretty on it's own and paste together 

p1<-ggplot(fmod17_predict, aes(x, predicted)) +
  geom_line() +
  geom_ribbon(data=fmod17_predict,aes(ymin = conf.low, ymax = conf.high), alpha = .1)+xlab("Provisioning Frequency (standardized)") + ylab("Predicted Foraging Time")+ggtitle("flock_size + provisioning + (1 | site)")+theme_classic()

p2<-ggplot(fmod17_predict3, aes(x, predicted)) +
  geom_line() +
  geom_ribbon(data=fmod17_predict3,aes(ymin = conf.low, ymax = conf.high), alpha = .1)+xlab("Flock Size (standardized)") + ylab(NULL)+theme_classic()

p12<-grid.arrange(p1,p2,ncol=2)

ggsave("outputs/foragingmodeloutput.png",plot=p12)
## Saving 7 x 5 in image
#make a plot that includes raw data and shows model fit:
#data = ForagePreen_wide
#this is Figure 3 in Manuscript

p1_points<-ggplot(fmod17_predict, aes(x, predicted)) +
  geom_line() +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .2,fill='#440154')+xlab("Provisioning frequency (standardized)") + ylab("Foraging time (seconds)")+ggtitle("flock_size + provisioning + (1 | site)")+theme_classic()+geom_point(data=ForagePreen_wide, aes(x=feedingperhour_s, y=Foraging),color='#ABABAB')

p2_points<-ggplot(fmod17_predict3, aes(x, predicted)) +
  geom_line() +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .2,fill='#440154')+xlab("Flock size (standardized)") + ylab("Foraging time (seconds)")+theme_classic()+geom_point(data=ForagePreen_wide, aes(x=flocksize_nearest_s, y=Foraging),color='#ABABAB')

fig3<-plot_grid(p1_points,p2_points, labels = c('A','B'))

#p12_points<-grid.arrange(p1_points,p2_points,ncol=2)
ggsave2("outputs/foragingmodeloutput_points.png",plot=fig3)
## Saving 7 x 5 in image

Make the same plots but without standardizing the predictor variables.

fmod17_unscaled <-glmmTMB(Foraging~flocksize_nearest+feedingperhour+(1|site)+offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=ForagePreen_wide)
summary(fmod17)
##  Family: nbinom2  ( log )
## Formula:          
## Foraging ~ flocksize_nearest_s + feedingperhour_s + (1 | site) +  
##     offset(log(duration))
## Zero inflation:            ~1
## Data: ForagePreen_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1567.2   1585.3   -777.6   1555.2      144 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev. 
##  site   (Intercept) 2.283e-09 4.779e-05
## Number of obs: 150, groups:  site, 5
## 
## Dispersion parameter for nbinom2 family (): 1.06 
## 
## Conditional model:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -0.97975    0.09296 -10.539   <2e-16 ***
## flocksize_nearest_s -0.24333    0.10301  -2.362   0.0182 *  
## feedingperhour_s    -0.18472    0.08670  -2.131   0.0331 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.0330     0.1878  -5.499 3.82e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(fmod17_unscaled) #p values look the same
##  Family: nbinom2  ( log )
## Formula:          Foraging ~ flocksize_nearest + feedingperhour + (1 | site) +  
##     offset(log(duration))
## Zero inflation:            ~1
## Data: ForagePreen_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1567.2   1585.3   -777.6   1555.2      144 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 2.52e-09 5.02e-05
## Number of obs: 150, groups:  site, 5
## 
## Dispersion parameter for nbinom2 family (): 1.06 
## 
## Conditional model:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -0.534394   0.160263  -3.334 0.000855 ***
## flocksize_nearest -0.004545   0.001918  -2.369 0.017822 *  
## feedingperhour    -0.344850   0.161854  -2.131 0.033119 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.0330     0.1878  -5.499 3.82e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fmod17_predict_unscaled <- ggpredict(fmod17_unscaled, terms = "feedingperhour")

fmod17_predict2_unscaled <- ggpredict(fmod17_unscaled, terms = "flocksize_nearest")

p1_points_unscaled<-ggplot(fmod17_predict_unscaled, aes(x, predicted)) +
  geom_line() +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .2,fill='#440154')+xlab("Provisioning Frequency") + ylab("Foraging Time (seconds)")+ggtitle("A).\n flock_size + provisioning + (1 | site)")+theme_classic()+geom_point(data=ForagePreen_wide, aes(x=feedingperhour, y=Foraging),color='#ABABAB')

p2_points_unscaled<-ggplot(fmod17_predict2_unscaled, aes(x, predicted)) +
  geom_line() +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .2,fill='#440154')+xlab("Flock Size") + ggtitle("B).")+ylab("Foraging Time (seconds)")+theme_classic()+geom_point(data=ForagePreen_wide, aes(x=flocksize_nearest, y=Foraging),color='#ABABAB')

p12_points_unscaled<-grid.arrange(p1_points_unscaled,p2_points_unscaled,ncol=2)

ggsave("foragingmodeloutput_points_unscaled.png",plot=p12_points_unscaled)
## Saving 7 x 5 in image
Diagnostics using Dharma

Testing the fit of fmod17 and fmod18- the two best fitting models that include site as a random effect.

First test the fit of fmod17

##  Family: nbinom2  ( log )
## Formula:          
## Foraging ~ flocksize_nearest_s + feedingperhour_s + (1 | site) +  
##     offset(log(duration))
## Zero inflation:            ~1
## Data: ForagePreen_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1567.2   1585.3   -777.6   1555.2      144 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev. 
##  site   (Intercept) 2.283e-09 4.779e-05
## Number of obs: 150, groups:  site, 5
## 
## Dispersion parameter for nbinom2 family (): 1.06 
## 
## Conditional model:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -0.97975    0.09296 -10.539   <2e-16 ***
## flocksize_nearest_s -0.24333    0.10301  -2.362   0.0182 *  
## feedingperhour_s    -0.18472    0.08670  -2.131   0.0331 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.0330     0.1878  -5.499 3.82e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.65756, p-value = 0.148
## alternative hypothesis: two.sided

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0019, p-value = 1
## alternative hypothesis: two.sided

Residuals look fine. Plotting residuals against the two predictors suggests that there might be some issues for one of them.

Now test fit for fmod18

##  Family: nbinom2  ( log )
## Formula:          
## Foraging ~ flocksize_nearest_s * feedingperhour_s + (1 | site) +  
##     offset(log(duration))
## Zero inflation:            ~1
## Data: ForagePreen_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1569.1   1590.2   -777.6   1555.1      143 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev. 
##  site   (Intercept) 2.117e-09 4.601e-05
## Number of obs: 150, groups:  site, 5
## 
## Dispersion parameter for nbinom2 family (): 1.07 
## 
## Conditional model:
##                                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                          -0.97710    0.09345 -10.455   <2e-16 ***
## flocksize_nearest_s                  -0.24949    0.10501  -2.376   0.0175 *  
## feedingperhour_s                     -0.20032    0.09952  -2.013   0.0441 *  
## flocksize_nearest_s:feedingperhour_s -0.05511    0.17401  -0.317   0.7515    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.0333     0.1879  -5.499 3.81e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.66608, p-value = 0.182
## alternative hypothesis: two.sided

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0021, p-value = 1
## alternative hypothesis: two.sided

The residuals look pretty okay for this model too but there are more problems when I plot against specific predictors.

Tables for Publication

Cand.models <- list("null + (1|site)" =fmod0.2, 
                    "provisioning + (1|site)" = fmod15,
                    "flock size + (1|site)" = fmod16,
                    "provisioning + flock size + (1|site)" = fmod17,
                    "provisioning * flock size + (1|site)" = fmod18)
selectionTable <- aictab(cand.set = Cand.models, second.ord = FALSE) #uses AIC 
row.names(selectionTable) <- NULL

kable(selectionTable,col.names = c("Model","K","AIC"     ,  "$\\Delta$ AIC","ModelLik","AICWt", "LogLik","Cum.Wt"),caption = "<b>Foraging Candidate Models<b>") %>%
 kable_styling(latex_options = "striped") %>%
  kable_classic(  lightable_options = "striped",full_width = F, html_font = "Cambria")
Foraging Candidate Models
Model K AIC \(\Delta\) AIC ModelLik AICWt LogLik Cum.Wt
provisioning + flock size + (1|site) 6 1567.219 0.000000 1.0000000 0.4603299 -777.6095 0.4603299
provisioning * flock size + (1|site) 7 1569.121 1.901545 0.3864425 0.1778910 -777.5603 0.6382209
flock size + (1|site) 5 1569.594 2.374851 0.3050054 0.1404031 -779.7969 0.7786240
provisioning + (1|site) 5 1569.842 2.623076 0.2694053 0.1240153 -779.9211 0.9026393
null + (1|site) 4 1570.326 3.107042 0.2115019 0.0973607 -781.1630 1.0000000
sf17<-summary(fmod17)
sumfmod17<-as.data.frame(sf17$coefficients[1])

kable(sumfmod17,col.names = c("Estimate","Std. error","z value","p-value"),caption = "<b>fmod17 summary<b>") %>%
 kable_styling(latex_options = "striped") %>%
  kable_classic(  lightable_options = "striped",full_width = F, html_font = "Cambria")
fmod17 summary
Estimate Std. error z value p-value
(Intercept) -0.9797530 0.0929608 -10.539422 0.0000000
flocksize_nearest_s -0.2433306 0.1030108 -2.362186 0.0181675
feedingperhour_s -0.1847193 0.0866983 -2.130599 0.0331222
sf18<-summary(fmod18)
sumfmod18<-as.data.frame(sf18$coefficients[1])

kable(sumfmod18,col.names = c("Estimate","Std. error","z value","p-value"),caption = "<b>fmod18 summary<b>") %>%
 kable_styling(latex_options = "striped") %>%
  kable_classic(  lightable_options = "striped",full_width = F, html_font = "Cambria")
fmod18 summary
Estimate Std. error z value p-value
(Intercept) -0.9770989 0.0934536 -10.4554490 0.0000000
flocksize_nearest_s -0.2494926 0.1050081 -2.3759363 0.0175045
feedingperhour_s -0.2003245 0.0995204 -2.0128988 0.0441253
flocksize_nearest_s:feedingperhour_s -0.0551059 0.1740110 -0.3166803 0.7514862
#try to merge these two summary tables 
sumfmod17$Model <- "fmod17"
sumfmod18$Model <- "fmod18"
sumfmod17$Predictors <- rownames(sumfmod17) #convert rownames to  column
sumfmod18$Predictors <- rownames(sumfmod18)#convert rownames to  column

fmods <- rbind(sumfmod17, sumfmod18)
rownames(fmods)<-NULL #remove rownames

fmods <- fmods %>%
  dplyr::select(Predictors, everything())

fmods$Model<-NULL
fmods <- fmods %>% 
       rename("pval" = "cond.Pr...z..") #renaming column 

#fmods$pval<-
 # dplyr::select(fmods,pval)%>% mutate_if(is.numeric, ~ case_when(. <0.00000001 ~ round(., 30), TRUE ~ round(.,5)))
#fmods <- fmods %>% 
   #    rename("pval" = "cond.Pr...z..") #renaming column 

#fmods  <- fmods  %>%
 # mutate(pval = case_when(
  # pval < 0.05 ~ paste0(as.character(pval), "*"),
  #  TRUE ~ paste0(pval)
 # ))

kable(fmods,col.names = c("Response/Predictors","Estimate","Std. error","z-value","p-value"),digits = c(5,5,5,5,20)) %>%
 kable_styling(latex_options = "striped") %>%
  kable_classic(  lightable_options = "striped",full_width = F, html_font = "Cambria")%>%
  pack_rows("Foraging", 1, 3) %>%
  pack_rows("Foraging", 4, 7)
Response/Predictors Estimate Std. error z-value p-value
Foraging
(Intercept) -0.97975 0.09296 -10.53942 0.00000000
flocksize_nearest_s -0.24333 0.10301 -2.36219 0.01816752
feedingperhour_s -0.18472 0.08670 -2.13060 0.03312218
Foraging
(Intercept) -0.97710 0.09345 -10.45545 0.00000000
flocksize_nearest_s -0.24949 0.10501 -2.37594 0.01750448
feedingperhour_s -0.20032 0.09952 -2.01290 0.04412528
flocksize_nearest_s:feedingperhour_s -0.05511 0.17401 -0.31668 0.75148618

Grooming models:

Main Q:

Does flock size, provisioning, both, or neither influence foraging time?

Grooming includes preening or bathing.

Exploratory Visualizations

Model selection:

Count data: Does poisson or negative binomial fit better? Build global model (two predictors + interactions) and compare AIC.

pm1_preen = glmmTMB(Preening~flocksize_nearest_s*feedingperhour_s + offset(log(duration)), ForagePreen_wide, family=poisson)

nbm1_preen = glmmTMB(Preening~flocksize_nearest_s*feedingperhour_s  +offset(log(duration)), ForagePreen_wide, family=nbinom2)

zipm1_preen = glmmTMB(Preening~flocksize_nearest_s*feedingperhour_s  +offset(log(duration)), zi=~1, ForagePreen_wide, family=poisson)

zinbm1_preen = glmmTMB(Preening~flocksize_nearest_s*feedingperhour_s +offset(log(duration)), zi=~1, ForagePreen_wide, family=nbinom2)


AICtab(pm1_preen,nbm1_preen,zipm1_preen,zinbm1_preen) #zero inflated is best here, continue with this model
##              dAIC    df
## zinbm1_preen     0.0 6 
## nbm1_preen      31.7 5 
## zipm1_preen  14544.6 5 
## pm1_preen    24584.0 4

Zero inflated negative binomial model is the best one so continue with that model structure.

##No random effects: 

#intercept only
pmod0.3 <-glmmTMB(Preening~1+ offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=ForagePreen_wide)
summary(pmod0.3)
##  Family: nbinom2  ( log )
## Formula:          Preening ~ 1 + offset(log(duration))
## Zero inflation:            ~1
## Data: ForagePreen_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1545.9   1554.9   -769.9   1539.9      147 
## 
## 
## Dispersion parameter for nbinom2 family (): 0.763 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -1.158      0.108  -10.72   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.1912     0.2037  -5.849 4.95e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#provisioning
pmod11 <-glmmTMB(Preening~feedingperhour_s+offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=ForagePreen_wide)
summary(pmod11)
##  Family: nbinom2  ( log )
## Formula:          Preening ~ feedingperhour_s + offset(log(duration))
## Zero inflation:            ~1
## Data: ForagePreen_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1547.2   1559.3   -769.6   1539.2      146 
## 
## 
## Dispersion parameter for nbinom2 family (): 0.766 
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -1.16616    0.10803 -10.795   <2e-16 ***
## feedingperhour_s  0.09753    0.12266   0.795    0.427    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.1919     0.2038  -5.848 4.98e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#flock size
pmod12 <-glmmTMB(Preening~flocksize_nearest_s+offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=ForagePreen_wide)
summary(pmod12)
##  Family: nbinom2  ( log )
## Formula:          Preening ~ flocksize_nearest_s + offset(log(duration))
## Zero inflation:            ~1
## Data: ForagePreen_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1545.0   1557.0   -768.5   1537.0      146 
## 
## 
## Dispersion parameter for nbinom2 family (): 0.783 
## 
## Conditional model:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -1.1765     0.1067 -11.022   <2e-16 ***
## flocksize_nearest_s   0.1638     0.1016   1.613    0.107    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.1870     0.2027  -5.856 4.74e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#provisioning + flock size
pmod13 <-glmmTMB(Preening~flocksize_nearest_s+feedingperhour_s+offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=ForagePreen_wide)
summary(pmod13)
##  Family: nbinom2  ( log )
## Formula:          
## Preening ~ flocksize_nearest_s + feedingperhour_s + offset(log(duration))
## Zero inflation:            ~1
## Data: ForagePreen_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1546.6   1561.6   -768.3   1536.6      145 
## 
## 
## Dispersion parameter for nbinom2 family (): 0.784 
## 
## Conditional model:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -1.18239    0.10677 -11.074   <2e-16 ***
## flocksize_nearest_s  0.15642    0.10063   1.554    0.120    
## feedingperhour_s     0.07797    0.11909   0.655    0.513    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.1875     0.2028  -5.855 4.76e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#provisioning x flock size 
pmod14 <-glmmTMB(Preening~flocksize_nearest_s*feedingperhour_s+offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=ForagePreen_wide)
summary(pmod14)
##  Family: nbinom2  ( log )
## Formula:          
## Preening ~ flocksize_nearest_s * feedingperhour_s + offset(log(duration))
## Zero inflation:            ~1
## Data: ForagePreen_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1548.6   1566.6   -768.3   1536.6      144 
## 
## 
## Dispersion parameter for nbinom2 family (): 0.784 
## 
## Conditional model:
##                                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                          -1.18402    0.10912 -10.851   <2e-16 ***
## flocksize_nearest_s                   0.15785    0.10248   1.540    0.123    
## feedingperhour_s                      0.08129    0.12786   0.636    0.525    
## flocksize_nearest_s:feedingperhour_s  0.01284    0.17902   0.072    0.943    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.1877     0.2029  -5.855 4.77e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AICtab(pmod0.3,pmod11,pmod12,pmod13,pmod14)
##         dAIC df
## pmod12  0.0  4 
## pmod0.3 0.9  3 
## pmod13  1.6  5 
## pmod11  2.2  4 
## pmod14  3.6  6
##Same models but with site as random effect 
#intercept only
pmod0.2 <-glmmTMB(Preening~1+ (1|site) + offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=ForagePreen_wide)
summary(pmod0.2)
##  Family: nbinom2  ( log )
## Formula:          Preening ~ 1 + (1 | site) + offset(log(duration))
## Zero inflation:            ~1
## Data: ForagePreen_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1545.8   1557.8   -768.9   1537.8      146 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.07115  0.2667  
## Number of obs: 150, groups:  site, 5
## 
## Dispersion parameter for nbinom2 family (): 0.802 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.2169     0.1624  -7.491 6.84e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.1857     0.2025  -5.856 4.73e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#provisioning
pmod15 <-glmmTMB(Preening~feedingperhour_s+(1|site)+offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=ForagePreen_wide)
summary(pmod15)
##  Family: nbinom2  ( log )
## Formula:          
## Preening ~ feedingperhour_s + (1 | site) + offset(log(duration))
## Zero inflation:            ~1
## Data: ForagePreen_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1547.7   1562.8   -768.9   1537.7      145 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.08273  0.2876  
## Number of obs: 150, groups:  site, 5
## 
## Dispersion parameter for nbinom2 family (): 0.805 
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -1.22063    0.17031  -7.167 7.65e-13 ***
## feedingperhour_s -0.03225    0.17501  -0.184    0.854    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.1851     0.2023  -5.857 4.72e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#flock size
pmod16 <-glmmTMB(Preening~flocksize_nearest_s+(1|site)+offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=ForagePreen_wide)
summary(pmod16)
##  Family: nbinom2  ( log )
## Formula:          
## Preening ~ flocksize_nearest_s + (1 | site) + offset(log(duration))
## Zero inflation:            ~1
## Data: ForagePreen_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1546.6   1561.7   -768.3   1536.6      145 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.03272  0.1809  
## Number of obs: 150, groups:  site, 5
## 
## Dispersion parameter for nbinom2 family (): 0.799 
## 
## Conditional model:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -1.2026     0.1384  -8.686   <2e-16 ***
## flocksize_nearest_s   0.1329     0.1219   1.091    0.275    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.1848     0.2022  -5.858 4.68e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#provisioning + flock size
pmod17 <-glmmTMB(Preening~flocksize_nearest_s+feedingperhour_s+(1|site)+offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=ForagePreen_wide)
summary(pmod17)
##  Family: nbinom2  ( log )
## Formula:          
## Preening ~ flocksize_nearest_s + feedingperhour_s + (1 | site) +  
##     offset(log(duration))
## Zero inflation:            ~1
## Data: ForagePreen_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1548.5   1566.6   -768.2   1536.5      144 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.01502  0.1226  
## Number of obs: 150, groups:  site, 5
## 
## Dispersion parameter for nbinom2 family (): 0.793 
## 
## Conditional model:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -1.19393    0.12748  -9.366   <2e-16 ***
## flocksize_nearest_s  0.14585    0.11664   1.250    0.211    
## feedingperhour_s     0.05801    0.15247   0.380    0.704    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.1861     0.2026  -5.855 4.77e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#provisioning x flock size 
pmod18 <-glmmTMB(Preening~flocksize_nearest_s*feedingperhour_s+(1|site)+offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=ForagePreen_wide)
summary(pmod18)
##  Family: nbinom2  ( log )
## Formula:          
## Preening ~ flocksize_nearest_s * feedingperhour_s + (1 | site) +  
##     offset(log(duration))
## Zero inflation:            ~1
## Data: ForagePreen_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1550.5   1571.6   -768.2   1536.5      143 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.01476  0.1215  
## Number of obs: 150, groups:  site, 5
## 
## Dispersion parameter for nbinom2 family (): 0.792 
## 
## Conditional model:
##                                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                          -1.194963   0.128629  -9.290   <2e-16 ***
## flocksize_nearest_s                   0.147157   0.118857   1.238    0.216    
## feedingperhour_s                      0.060848   0.161153   0.378    0.706    
## flocksize_nearest_s:feedingperhour_s  0.009679   0.179544   0.054    0.957    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.1862     0.2026  -5.855 4.78e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AICtab(pmod0.2,pmod15,pmod16,pmod17,pmod18)
##         dAIC df
## pmod0.2 0.0  4 
## pmod16  0.8  5 
## pmod15  2.0  5 
## pmod17  2.7  6 
## pmod18  4.7  7

The top performing models are within 2 AIC values of the intercept only null model.

Model diagnostics

Not needed because the intercept only model is the best fitting model?

Tables for Publication

Cand.models_groom <- list("null + (1|site)" =pmod0.2, 
                    "provisioning + (1|site)" = pmod15,
                    "flock size + (1|site)" = pmod16,
                    "provisioning + flock size + (1|site)" = pmod17,
                    "provisioning * flock size + (1|site)" = pmod18)
selectionTable_groom <- aictab(cand.set = Cand.models_groom, second.ord = FALSE) #uses AIC 
row.names(selectionTable_groom) <- NULL

kable(selectionTable_groom,col.names = c("Model","K","AIC"     ,  "$\\Delta$ AIC","ModelLik","AICWt", "LogLik","Cum.Wt"),caption = "<b>Grooming Candidate Models<b>") %>%
 kable_styling(latex_options = "striped") %>%
  kable_classic(  lightable_options = "striped",full_width = F, html_font = "Cambria")
Grooming Candidate Models
Model K AIC \(\Delta\) AIC ModelLik AICWt LogLik Cum.Wt
null + (1|site) 4 1545.776 0.0000000 1.0000000 0.4200894 -768.8879 0.4200894
flock size + (1|site) 5 1546.623 0.8474312 0.6546100 0.2749948 -768.3116 0.6950842
provisioning + (1|site) 5 1547.742 1.9657806 0.3742279 0.1572092 -768.8708 0.8522934
provisioning + flock size + (1|site) 6 1548.494 2.7177867 0.2569450 0.1079399 -768.2468 0.9602332
provisioning * flock size + (1|site) 7 1550.491 4.7148714 0.0946627 0.0397668 -768.2454 1.0000000

Vigilance models:

Main Q:

Does flock size, provisioning, both, or neither influence vigilance time?

Exploratory Visualizations

Model selection:

Count data: Does poisson or negative binomial fit better? Build global model (two predictors + interactions) and compare AIC.

pm1_vig = glmmTMB(Vigilant~flocksize_nearest_s*feedingperhour_s + offset(log(duration)), Vig_wide, family=poisson)

nbm1_vig = glmmTMB(Vigilant~flocksize_nearest_s*feedingperhour_s  +offset(log(duration)), Vig_wide, family=nbinom2)

zipm1_vig = glmmTMB(Vigilant~flocksize_nearest_s*feedingperhour_s  +offset(log(duration)), zi=~1, Vig_wide, family=poisson)

zinbm1_vig = glmmTMB(Vigilant~flocksize_nearest_s*feedingperhour_s +offset(log(duration)), zi=~1, Vig_wide, family=nbinom2)


AICtab(pm1_vig,nbm1_vig,zipm1_vig,zinbm1_vig) #zero inflated is best here, continue with this model
##            dAIC    df
## zinbm1_vig     0.0 6 
## nbm1_vig      14.0 5 
## zipm1_vig  17292.7 5 
## pm1_vig    19778.4 4

Zero inflated negative binomial model is the best one so continue with that model structure.

##No random effects: 

#intercept only
vmod0.3 <-glmmTMB(Vigilant~1+ offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=Vig_wide)
summary(vmod0.3)
##  Family: nbinom2  ( log )
## Formula:          Vigilant ~ 1 + offset(log(duration))
## Zero inflation:            ~1
## Data: Vig_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1795.1   1804.1   -894.5   1789.1      147 
## 
## 
## Dispersion parameter for nbinom2 family (): 1.04 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.1186     0.0831  -13.46   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.7527     0.3674  -7.493 6.72e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#provisioning
vmod11 <-glmmTMB(Vigilant~feedingperhour_s+offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=Vig_wide)
summary(vmod11)
##  Family: nbinom2  ( log )
## Formula:          Vigilant ~ feedingperhour_s + offset(log(duration))
## Zero inflation:            ~1
## Data: Vig_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1795.2   1807.3   -893.6   1787.2      146 
## 
## 
## Dispersion parameter for nbinom2 family (): 1.05 
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -1.12883    0.08280 -13.634   <2e-16 ***
## feedingperhour_s  0.11559    0.08451   1.368    0.171    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.7571     0.3691  -7.471 7.98e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#flock size
vmod12 <-glmmTMB(Vigilant~flocksize_nearest_s+offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=Vig_wide)
summary(vmod12)
##  Family: nbinom2  ( log )
## Formula:          Vigilant ~ flocksize_nearest_s + offset(log(duration))
## Zero inflation:            ~1
## Data: Vig_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1796.7   1808.7   -894.3   1788.7      146 
## 
## 
## Dispersion parameter for nbinom2 family (): 1.05 
## 
## Conditional model:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -1.11978    0.08298 -13.494   <2e-16 ***
## flocksize_nearest_s -0.05246    0.08319  -0.631    0.528    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.7504     0.3665  -7.504 6.21e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#provisioning + flock size
vmod13 <-glmmTMB(Vigilant~flocksize_nearest_s+feedingperhour_s+offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=Vig_wide)
summary(vmod13)
##  Family: nbinom2  ( log )
## Formula:          
## Vigilant ~ flocksize_nearest_s + feedingperhour_s + offset(log(duration))
## Zero inflation:            ~1
## Data: Vig_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1796.6   1811.7   -893.3   1786.6      145 
## 
## 
## Dispersion parameter for nbinom2 family (): 1.06 
## 
## Conditional model:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -1.13097    0.08263 -13.688   <2e-16 ***
## flocksize_nearest_s -0.06498    0.08119  -0.800    0.423    
## feedingperhour_s     0.12392    0.08560   1.448    0.148    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.7546     0.3682  -7.482 7.33e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#provisioning x flock size 
vmod14 <-glmmTMB(Vigilant~flocksize_nearest_s*feedingperhour_s+offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=Vig_wide)
summary(vmod14)
##  Family: nbinom2  ( log )
## Formula:          
## Vigilant ~ flocksize_nearest_s * feedingperhour_s + offset(log(duration))
## Zero inflation:            ~1
## Data: Vig_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1798.5   1816.6   -893.2   1786.5      144 
## 
## 
## Dispersion parameter for nbinom2 family (): 1.06 
## 
## Conditional model:
##                                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                          -1.13733    0.08423 -13.503   <2e-16 ***
## flocksize_nearest_s                  -0.05786    0.08320  -0.695    0.487    
## feedingperhour_s                      0.13835    0.09467   1.461    0.144    
## flocksize_nearest_s:feedingperhour_s  0.05327    0.14684   0.363    0.717    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.7558     0.3687  -7.474 7.75e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AICtab(vmod0.3,vmod11,vmod12,vmod13,vmod14)
##         dAIC df
## vmod0.3 0.0  3 
## vmod11  0.2  4 
## vmod13  1.5  5 
## vmod12  1.6  4 
## vmod14  3.4  6
##Same models but with site as random effect 
#intercept only
vmod0.2 <-glmmTMB(Vigilant~1+ (1|site) + offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=Vig_wide)
summary(vmod0.2)
##  Family: nbinom2  ( log )
## Formula:          Vigilant ~ 1 + (1 | site) + offset(log(duration))
## Zero inflation:            ~1
## Data: Vig_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1796.7   1808.8   -894.4   1788.7      146 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.01613  0.127   
## Number of obs: 150, groups:  site, 5
## 
## Dispersion parameter for nbinom2 family (): 1.06 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.1253     0.1008  -11.16   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.7493     0.3662  -7.508 5.99e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#provisioning
vmod15 <-glmmTMB(Vigilant~feedingperhour_s+(1|site)+offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=Vig_wide)
summary(vmod15)
##  Family: nbinom2  ( log )
## Formula:          
## Vigilant ~ feedingperhour_s + (1 | site) + offset(log(duration))
## Zero inflation:            ~1
## Data: Vig_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1796.9   1812.0   -893.5   1786.9      145 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.0167   0.1292  
## Number of obs: 150, groups:  site, 5
## 
## Dispersion parameter for nbinom2 family (): 1.07 
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -1.13354    0.10079 -11.247   <2e-16 ***
## feedingperhour_s  0.12838    0.09711   1.322    0.186    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.7515     0.3671  -7.495 6.61e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#flock size
vmod16 <-glmmTMB(Vigilant~flocksize_nearest_s+(1|site)+offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=Vig_wide)
summary(vmod16)
##  Family: nbinom2  ( log )
## Formula:          
## Vigilant ~ flocksize_nearest_s + (1 | site) + offset(log(duration))
## Zero inflation:            ~1
## Data: Vig_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1798.4   1813.4   -894.2   1788.4      145 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.01607  0.1268  
## Number of obs: 150, groups:  site, 5
## 
## Dispersion parameter for nbinom2 family (): 1.06 
## 
## Conditional model:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -1.12738    0.10085  -11.18   <2e-16 ***
## flocksize_nearest_s -0.05547    0.09087   -0.61    0.542    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.7480     0.3657  -7.514 5.73e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#provisioning + flock size
vmod17 <-glmmTMB(Vigilant~flocksize_nearest_s+feedingperhour_s+(1|site)+offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=Vig_wide)
summary(vmod17)
##  Family: nbinom2  ( log )
## Formula:          
## Vigilant ~ flocksize_nearest_s + feedingperhour_s + (1 | site) +  
##     offset(log(duration))
## Zero inflation:            ~1
## Data: Vig_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1798.5   1816.6   -893.2   1786.5      144 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.01066  0.1032  
## Number of obs: 150, groups:  site, 5
## 
## Dispersion parameter for nbinom2 family (): 1.07 
## 
## Conditional model:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -1.13444    0.09477 -11.971   <2e-16 ***
## flocksize_nearest_s -0.05885    0.08814  -0.668    0.504    
## feedingperhour_s     0.12762    0.09276   1.376    0.169    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.7515     0.3671  -7.495 6.64e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#provisioning x flock size 
vmod18 <-glmmTMB(Vigilant~flocksize_nearest_s*feedingperhour_s+(1|site)+offset(log(duration)),
                    family=nbinom2,
                    ziformula=~1,
                    data=Vig_wide)
summary(vmod18)
##  Family: nbinom2  ( log )
## Formula:          
## Vigilant ~ flocksize_nearest_s * feedingperhour_s + (1 | site) +  
##     offset(log(duration))
## Zero inflation:            ~1
## Data: Vig_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##   1800.4   1821.5   -893.2   1786.4      143 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.007085 0.08417 
## Number of obs: 150, groups:  site, 5
## 
## Dispersion parameter for nbinom2 family (): 1.06 
## 
## Conditional model:
##                                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                          -1.13790    0.09207 -12.359   <2e-16 ***
## flocksize_nearest_s                  -0.05544    0.08769  -0.632    0.527    
## feedingperhour_s                      0.13647    0.09856   1.385    0.166    
## flocksize_nearest_s:feedingperhour_s  0.03804    0.16156   0.235    0.814    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.7533     0.3679  -7.483 7.24e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AICtab(vmod0.2,vmod15,vmod16,vmod17,vmod18)
##         dAIC df
## vmod0.2 0.0  4 
## vmod15  0.2  5 
## vmod16  1.6  5 
## vmod17  1.8  6 
## vmod18  3.7  7

The top performing models are within 2 AIC values of the intercept only null model.

Model diagnostics

Not needed because the intercept only model is the best fitting model?

Tables for Publication

Cand.models_vig <- list("null + (1|site)" =vmod0.2, 
                    "provisioning + (1|site)" = vmod15,
                    "flock size + (1|site)" = vmod16,
                    "provisioning + flock size + (1|site)" = vmod17,
                    "provisioning * flock size + (1|site)" = vmod18)
selectionTable_vig <- aictab(cand.set = Cand.models_vig, second.ord = FALSE) #uses AIC 
row.names(selectionTable_vig) <- NULL

kable(selectionTable_vig,col.names = c("Model","K","AIC"     ,  "$\\Delta$ AIC","ModelLik","AICWt", "LogLik","Cum.Wt"),caption = "<b>Vigilance Candidate Models<b>") %>%
 kable_styling(latex_options = "striped") %>%
  kable_classic(  lightable_options = "striped",full_width = F, html_font = "Cambria")
Vigilance Candidate Models
Model K AIC \(\Delta\) AIC ModelLik AICWt LogLik Cum.Wt
null + (1|site) 4 1796.749 0.0000000 1.0000000 0.3405399 -894.3743 0.3405399
provisioning + (1|site) 5 1796.914 0.1650839 0.9207728 0.3135599 -893.4569 0.6540997
flock size + (1|site) 5 1798.384 1.6356869 0.4413825 0.1503083 -894.1922 0.8044081
provisioning + flock size + (1|site) 6 1798.500 1.7508891 0.4166767 0.1418950 -893.2498 0.9463031
provisioning * flock size + (1|site) 7 1800.443 3.6943541 0.1576817 0.0536969 -893.2215 1.0000000